#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<omp.h>
#include<time.h>
#include<ctype.h>

void comprovacao(int J_ORDER, double x[J_ORDER], int J_ROW_TEST, double **matriz, double b[J_ORDER], int it) {
	int i;
	double valor = 0;
        
    for(i = 0; i < J_ORDER; i++) {
    	valor += matriz[J_ROW_TEST][i]*x[i];
    }
        
    printf("---------------------------------------\n");
    printf("Iterations: %d\n", it);
    printf("RowTest: %d => [%lf] =? %lf\n", J_ROW_TEST, valor, b[J_ROW_TEST]);
    printf("---------------------------------------\n");
        
}

int verifica(int J_ORDER, double **aEstrela) {
    int i, j;
    double x[J_ORDER], max;

    for(i = 0; i < J_ORDER; i++) {
        x[i] = 0;
    }

    for(i = 0; i < J_ORDER; i++) {
        for(j = 0; j < J_ORDER; j++) {
            if(i != j){
                x[i] += aEstrela[i][j];
            }
        }
    }
    max = x[0];

    for(i = 1; i < J_ORDER; i++) {
       if(max < x[i])
           max = x[i];
    }

    if(max < 1)
        return 1;

    else
        return 0;
}

int main(int argc, char *argv[]) {
	
	int J_ORDER, J_ROW_TEST, J_ITE_MAX, i, j;
	double J_ERROR, **matriz, **aEstrela;
	char narq[15];
	clock_t inicio = 0, fim = 0;
    inicio = clock();          
  	FILE *pFile;

  	printf("Nome do arquivo: ");
  	scanf("%s", narq);

  	pFile = fopen ( narq , "rb" );
  	if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

  	rewind (pFile); //volta para o início do arquivo
  	
  	fscanf(pFile, "%d", &J_ORDER);
  	fscanf(pFile, "%d", &J_ROW_TEST);
  	fscanf(pFile, "%lf", &J_ERROR);
  	fscanf(pFile, "%d", &J_ITE_MAX);

	double b[J_ORDER], x0[J_ORDER];
  
  	matriz = malloc( J_ORDER * sizeof(double*) );
	#pragma omp parallel for
  	for(i = 0; i < J_ORDER; i++) {
    	matriz[i] = malloc( J_ORDER * sizeof(double) );
  	}
  	
  	aEstrela = malloc( J_ORDER * sizeof(double*) );
  	#pragma omp parallel for
    for(i = 0; i < J_ORDER; i++) {
        aEstrela[i] = malloc( J_ORDER * sizeof(double) );
    }	
  	
  	for(i = 0; i < J_ORDER; i++) {
    	for(j = 0; j < J_ORDER; j++) {
        	fscanf(pFile, "%lf", &matriz[i][j]);
        }
  	}
  	
  	for(i = 0; i < J_ORDER; i++) {
        fscanf(pFile, "%lf", &b[i]);
        x0[i] = 0;
  	}
  	
  	for(i = 0; i < J_ORDER; i++) {
    	for(j = 0; j < J_ORDER; j++) {
        	if(i != j) {
            	aEstrela[i][j] = (matriz[i][j] / matriz[i][i])*(-1);
        	}
        }
    }
    
    if(verifica(J_ORDER, aEstrela) != 1) {
        return;
    }
    else {
    	int it = 0;//numero de iterações
        double bEstrela[J_ORDER], **x;
        
        x = malloc( J_ITE_MAX * sizeof(double*) );
        for(i = 0; i < J_ITE_MAX; i++) {
                x[i] = malloc( J_ORDER * sizeof(double) );
        }
        
        for(i = 0; i < J_ORDER; i++){
        	x[0][i] = x0[i];
            bEstrela[i] = b[i]/matriz[i][i];
        }
    
		double e = 0, max1, max2;// e = erro relativo\/max1 auxiliar para ver qual a maior diferença entre Xk+1 e Xk e max2 é o maior X

		do {
			it++;
			//Formula iterativa;
			int contx = 0;
			
			#pragma omp parallel for firstprivate(j,x)
			for(i = 0; i < J_ORDER; i++) {  
		    	for(j = 0; j < J_ORDER; j++) {
		        	if(i != j) {
		            	x[it][i] += (aEstrela[i][j] * x[it-1][j]);
					}	 
				}
		        x[it][i] = x[it][i] + bEstrela[i];
		    }

		    //acha o maior valor de Xk+1 - Xk (1 por 0)
		    max1 = fabs(x[it][1] - x[it-1][1]);
		    max2 = fabs(x[it][1]);

		    for(i = 1; i < J_ORDER; i++) {
		    	if(max1 < fabs (x[it][i] - x[it-1][i]) ) {
		        	max1 = fabs (x[it][i] - x[it-1][i]);
		    	}

		    	if(max2 < x[it][i] ) {
		        	max2 = fabs(x[it][i]);
		    	}

		    }

		    e = max1/max2;

		}while(e > J_ERROR);
	
		comprovacao(J_ORDER, x[it], J_ROW_TEST, matriz, b, it);
		
		#pragma omp parallel for 
        for(i = 0; i < J_ITE_MAX; i++) {
        	free(x[i]);
        }
        free(x);
	}
	
	#pragma omp parallel for 
    for(i = 0; i < J_ORDER; i++) {
     	free(aEstrela[i]);
    }
    free(aEstrela);
	
	#pragma omp parallel for 
    for(i = 0; i < J_ORDER; i++) {
    	free(matriz[i]);
    }
    free(matriz);
    
 	fim  = clock();
	double  tempo = (fim - inicio)/(CLOCKS_PER_SEC/1000);
	printf("tempo = %g\n", tempo);

	fclose (pFile);

  	return 0;
	
}
